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ABSTRACT 


A nonlinear simulation is developed to model the longitudinal motion of a vehicle in 
hypersonic flight The equations of motion pertinent to this study are presented. Analytic 
expressions for the aerodynamic forces acting on a hypersonic vehicle which were obtained 
from Newtonian Impact Theory are further developed. The control surface forces are 
further examined to incorporate vehicle elastic motion. The purpose is to establish feasible 
equations of motion which combine rigid body, elastic and aeropropulsive dynamics for 
use in nonlinear simulations. The software package SIMULINK® is used to implement the 
simulation. Also discussed are issues needing additional attention and potential problems 
associated with the implementation (with proposed solutions). 
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1. INTRODUCTION 


The object of the Hypersonic Vehicle Simulation ASUHS1 is to integrate the equations of 
motion for a flexible hypersonic vehicle in atmospheric flight. The derivation of the 
equations of motion is fully documented in Reference [1]. An aeroelastic model is 
presented in Reference [2]. This class of vehicle has been shown to exhibit a high degree 
of coupling among the aerodynamic, propulsive, and elastic effects. It also exhibits highly 
unstable pitch behavior [4]. 

What follows is a discussion of ASUHS1, which simulates the dynamics of the hypersonic 
vehicle model using the simulation software package SIMULINK® [3]. The nonlinear 
simulation developed in this report will serve as a tool to better understand hypersonic 
vehicle dynamics. For this simulation, only motion in a vertical plane (symmetric flight) is 
considered. 

Reference [2] includes analytical expressions for the linearized aerodynamic and propulsive 
forces on the vehicle. The total aerodynamic forces on the lower forebody, however, are 
expressed in terms of (nonlinear) integral equations. In this report, these integral equations 
are solved in closed form and used in the simulation. In addition, the equations for the 
aerodynamic control surface forces and moments are modified to incorporate the effect of 
vehicle bending not included in [2]. 

The organization of this report is as follows: Chapter 2 first outlines the structure of 
ASUHS1. Chapter 3 presents the equations of motion for this class of vehicle and 
simplifies them for motion restricted to the vertical plane. In Chapter 4 the aeropropulsive 
force equations are presented. The computational procedure is discussed in Chapter 5, 
with numerical data given in Chapter 6. Potential problems of which the user should be 
aware complete this report. 


2. STRUCTURE OF SIMULATION 


The simulation structure is best explained by referring to Figure 1. The two basic tasks of 
ASUHS 1 are to compute the aeropropulsive forces acting on the vehicle and to integrate the 
equations of motion for this class of vehicle. The first of the above tasks is carried out in 
the block FORCES, while the second is done in the block EOM. Each of these blocks 
consists of several sub-blocks, all of which utilize common elements found in the 
SIMULINK® library and also coded functions (M-files) found in the local directory. 
These MATLAB® functions are fully documented in Appendix A. 

The initialization procedure for ASUHS 1 consists of running two script files, inpar and 
initstate. The first file initializes the parameter values associated with the spherical rotating 
earth and atmospheric models. These values will generally remain the same regardless of 
the study being conducted. The vehicle's length, elastic and mass properties, and engine 
parameters, which may vary from study to study, are also defined in this script file. The 
second file initializes the initial conditions and control settings. These values may be 
altered as well, according to the desired flight condition. 

The aeropropulsive forces (and moments) are categorized as: l)Aerodynamic forces on the 
vehicle’s lower forebody, 2)Aerodynamic forces on the vehicle’s control surface, 3)Engine 
thrust forces, 4)Extemal nozzle exhaust forces. These forces are computed from separate 
functions of vehicle states and controls. Calculation of the these forces includes the effect 
of the vehicle’s elastic motion, which changes the angle of the lower forebody and lower 
aft-body, as well as control surface position relative to the vehicle’s center of gravity. The 
forces are summed just before leaving the FORCES block. 

The block EOM makes use of the vehicle equations of motion over a spherical rotating 
earth [1]. This block utilizes the above forces and moments, plus the gravitational "force." 
The state derivatives are calculated and numerically integrated. The method of numerical 
integration is chosen by the user from SIMULINK® ’s simulation parameters. 
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Figure 1. Hypersonic Vehicle Simulation block structure. 








3. EQUATIONS OF MOTION 


An integrated model for elastic hypersonic vehicles was rigorously developed in Reference 
[1]. This approach accounted for the interaction between rigid body motion, elastic 
deformation, fluid flow, rotating machinery, wind and a spherical rotating earth. What 
resulted were 12+n equations of motion: the 3 force equations + 3 moment equations + n 
elastic deformation equations + 3 kinematic equations + 3 Euler angle equations. In this 
simulation only symmetric motion (i.e., motion in the vertical plane) is currently 
considered. 

For this initial development of the simulation, many of the interactions accounted for in [ 1 ] 
are neglected. For example, the effects of rotating machinery are neglected. Forces 
associated with mass flow effects not directly contributing to engine thrust are, at this time, 
considered negligible in comparison to the aeropropulsive forces; the corresponding 
moments are also neglected. A no-wind assumption eliminates all the wind velocity 
components. Finally, the vehicle is assumed to possess a plane of symmetry that passes 
through the body x-axis and z-axis with constant moments of inertia. These assumptions 
simplify the force equations (4.1-17) to (4.1-19), and moment equations (4.2-8) to (4.2- 
10) in [1]. 

To further reduce the 12+n equations of motion, aeropropulsive roll and yaw moments and 
body y-axis force (lateral force) are assumed absent [2]. Secondly, in ASUHS1 a wings- 
level flight condition is currently assumed. These two facts eliminate roll motion (<()=p=0) 
and yaw rate (r=0) from the present simulation. In addition, only one elastic degree of 
freedom is considered. This leaves three force equations, one moment equation and one 
elastic deformation equation as the field equations to be implemented. They are the body x- 
axis force equation (4.1-17) [1], y-axis force equation (4.1-18) [1], z-axis force equation 
(4.1-19) [1], y-axis moment equation (4.2-9) [1] and the first generalized elastic 
deformation equation (4.3-2) [1], These equations are listed below: 


u = -wq - 0 ) e [w(T 2 iCosX. - T 2 3 sinA,) - v(T 3 icosA. - T 33 sin\)] 
- RC0?[T, 1 SilOX + T 13 C0S2X] + go(|f Ti3 + ^ 

v = -0)e[u(T 3 icosX - T 33 sinX.) - w(TncosX - T^sinX,)] 

- + T 23 cos2X] + go(^fT 2 3 + £ 


(3.1) 


(3.2) 
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where. 


w = uq - co e [v(TiicosA, - Ti 3 sinX) - u(T 2 iC 0 sX. - T 2 3 sinA,)] 
- Rm?[T 3 iSil^i + T 3 ! cos 2 aJ + g,^) 2 T 33 + 2 . 

lyy 


ii = -« 5 n - 2 CiC0iT| + 


Qei 

m El 


(3.3) 

(3.4) 

(3.5) 


C 0 e = Earth rotation rate about its axis 
R e = Earth radius 

R = Distance from vehicle c.g. to Earth's center (R=R e +h) 
g 0 = Gravitational acceleration at the Earth’s surface 
X = Latitude of vehicle position 
m = Vehicle mass (per unit width of the vehicle) 

lyy = Vehicle moment of inertia about body y-axis (per unit width of the 
vehicle) 

C0j = Undamped natural frequency of 1 st generalized elastic degree of freedom 
= Damping ratio of 1 st generalized elastic degree of freedom 

m El = Generalized mass of 1 st generalized elastic degree of freedom (per unit 
width of the vehicle) 

T = Co-ordinate transformation matrix from vehicle-carrying frame (X v , Y v , 
Z v ) to body frame (X B , Y B , Z B ), see Figure 2. The entries of the co- 
ordinate transformation matrix T may be written in terms of the quaternion 
components described in [1], or in terms of the roll, pitch and yaw angles 
(<]), 0, V|/). Both representations for T are shown below: 



(3?-pl-pi+pi 

2(PiP 2 +P 3 p4) 

2(piP 3 -p 2 P4) 

T = 

2(Pip 2 -p 3 P4) 

-p?+p2-p 3 +pj 

2(p 2 P 3 +Plp4) 


2(PiP3+P2P4) 

2(p2p3-PlP 4 ) 

-P?-P2+P 3 +P4 . 


COS0COS\y 

cos0sin\jt 

-sin0 

T = 

sin<t>sin0cosY 

sin())sin0sinv|/ 

sin<)»cos0 


- cos<|)sin\j/ 

+ COS<j)COS\|/ 


cos<|)sin0cos\|r 

cos())sin0sin\|/ 

COS(()COS0 


+ sin<t>sin\|/ 

- sin<))cos\|r 



The following equations may be used to determine the quaternions from the Euler angles: 

^ sin(6)sin(v|/) 

4 cos( 0 / 2 )cos( 0 / 2 ) 

sin(0)cos(\|//2) 

P2= 

2cos(0/2) 

_ sin(\j/)cos(0/2) 

P3 2cos(\j//2) 

P4 = cos(0/2)cos(\|//2) 


N ' 



Figure 2. Vehicle trajectory variables and coordinate frames [1]. 



The implication of symmetric flight constrains motion to a vertical plane fixed in an inertial 
reference (in this case, Earth fixed center). But we will require for now that the vehicle's 
lateral velocity in the "E" frame (see Figure 2) equal zero (v=0) for all time. Obviously, the 
initial condition and constraint equation which meet the requirements of symmetric flight 
are: 

v(t=0) = v 0 = 0 

v(t) = 0, V t>0 (3.6) 


To implement equation (3.6) requires that the applied lateral force, Y, cancel the Coriolis 
and centripetal acceleration terms in equation (3.2). (Note that the gravity term in equation 
(3.2) is zero because T23=0 for wings level.) In conventional aircraft dynamics these terms 
are usually neglected since the vehicle speed is relatively small. This is not necessarily the 
case in hypersonic flight; consequently the Coriolis and centripetal acceleration terms are 
taken into account. In light of this unorthodox circumstance, Y is utilized in this initial 
model as a constraining force to ensure symmetric flight, and the lateral force equation is 
eliminated. The simplified force equations are: 


u = -w(q + co e T 2 icosX.) - RtoflTii ^^ 



+ Ti3COS 2 A.] 


w = u(q + 0 ) e T 2 icosA.) - R(o|[T 3 i s * 1 ] 2 X + T 33 Cos 2 X] 



(3.7) 


(3.8) 


The aeropropulsive forces and moments (X, Z, M and Qei) which appear in equations 
(3.7), (3.8), (3.4) and (3.5), respectively, are expressed in terms of the dynamic variables 
and controls. These expressions are examined in Chapter 4. The twelve dynamic variables 
to be integrated with three control inputs are: 


h = Geometric altitude - [ft] 

u = Vehicle forward speed - [ft/s] (relative to "E" frame) 
w = Vehicle sink rate - [ft/s] (relative to "E" frame) 
q = Vehicle rigid body pitch rate - [rad/s] 

T| = First generalized elastic co-ordinate - [-] 

B = First generalized elastic co-ordinate rate - [ 1/s] 

X = Latitude of vehicle position - [rad] 


x = Longitude of vehicle position - [rad] 

Pi = First quaternion component - [-] 

P2 = Second quaternion component - [-] 

P3 = Third quaternion component - [-] 

P4 = Fourth quaternion component - [-] 

5 = Pitch control surface deflection - [rad] 

Aq = Diffuser area ratio control - [-] 

T 0 = Combustor total temperature change control - [R°] 

In addition, three kinematic equations, one differential identity and four quaternion 
equations are necessary to describe the rates of change of the altitude, latitude, longitude, 
first generalized elastic co-ordinate and quaternions: 


h = -[T13U + T33W] ( 3 . 9 ) 

A-=|L[Tiiu + T 3 iw] ( 3 . 10 ) 

t = S&&[t 12U + T 3 2w] ( 3 . 11 ) 

Iv 

(3.i2) 

Pi =^-(PvP 4 * q v p3 +r v p 2 ) (3.13) 

P2 = 2 (p v p 3 + q v p 4 - r vPi) (3.14) 

P3 = 2 ("P v p 2 ■*" QvPi + r vP 4 ) (3. 15) 

P 4 = y(-PvPi -qvP2-r v p3) (3.16) 


where, Pv = r [uT 12 Ti 3 tanA, + w(Ti 3 T 32 tanA. + Ti 2 T 3 i - T n T 32 )] 

+ 0) e (Ti3sinX - TucosX) 

q v = q + ^- [u(ThT 22 - Ti 2 T 2 j)+ w( T 22 T31 - T 2 iT3 2 )] - 0) e T 2 iCosX 

Fv = R t u ( T l2 T 33tanA. + T11T32 - T12T31) + \vT32T33tanA.] 

+ 0) e (T33sinA. - T31COSA.) 

The block EOM determines the right-hand sides of equations ( 3 . 4 )-( 3 . 16 ) for time 
integration. 
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4. AEROPROPULSIVE FORCES AND MOMENTS 
4.1. Introduction 

Analytical expressions for the aerodynamic and propulsive forces acting on a hypersonic 
vehicle have been developed in [2]. Also, a linearized model has been obtained for modal 
analysis and control law design [2, 4]. For implementation purposes, the integral 
expressions in the model may be solved in closed form. 

For completeness, the expressions for the aeropropulsive forces are rewritten here from 
[2]. Recall they are composed of four separate constituents: l)Aerodynamic forces on the 
forebody, 2)Aerodynamic forces on the control surface, 3)Engine thrust, 4)External nozzle 
exhaust. 


. f 

-pnl 

Jo 


1) X A =-{P«li+q~C P n I sin^L^Odsil-sinCxi+AxiTi) 


■i 


Z A =P 00 L-{PJi+q 00 C pn | sin 2 0Lf(si)dsi}-cos(ti+Ati'n) 


(4.1-1) 

(4.1-2) 


P 1^ f h 

MA={-y^+qo.C pn | sisin 2 0Lf< s i)dsi}-[(h-z C g)sin(Xi+Axi'n)+(Li-x C g)cos(Xi-i-AxiTi)]- 

^ Jo 


{ PJi+q„C pn | sin 2 0 L f<si)dsi ) - P 00 L(x cg -^) 


PJi / ‘ 1 ' 


Q A = Ax i • { — — i-q^Cpn I sisin 2 0 L K s i)dsi } 


■ 

Jo 


2) X cs =-q„C pn sin 2 (a+5)sin(8)^- 


Zcs=-q~C pn sin 2 (a+5)cos(8)^ s - 


M cs — X cs z cs -Z cs x c 


(4.1-3) 

(4.1-4) 

(4.1-5) 

(4.1-6) 

(4.1-7) 


Qcs=0 


(4.1-8) 



(4.1-9) 


3) X T ={[YP e M|+(P e -PJ]- 


[TfPiMfoPt-PJ] Ae 
AqAn b 


Z T = 0 (4.1-10) 

Mj =Xx(h-Zc g ) (4. 1 - 1 1 ) 


Qr=0 


(4.1-12) 


4) X E =PJd 


Pln(p) 


(p-i) 


sin(x 2 +Ax 2 Ti) 


Ze= .p .ME 


cos(x 2 +Ax 2 ti) 


Me =P„ !bri^ (P| 


(P-1) 




, MS 

4 / — \ 


Ml M. 


Qe — A'^P.Je”,"^” 


(p-1) 


fi-MEl 


(p-l)J 


where, P = ^ 


r 2 = (h-Zcg)sin(x 2 +Ax 2 T|) - (Li-x cg )cos(x 2 +Ax 2 ri) 


(4.1-13) 

(4.1-14) 

(4.1-15) 

(4.1-16) 


4.2. Algebraic Solutions 

4.2.1 Explicit Expressions for Local Flow Deflection Angles 

The pressure distribution on a hypersonic vehicle is strongly dependent upon the air flow 
deflection angle [2], (See Figure 3.) The lower forebody of the vehicle deflects the air 
flow an angle 0. Due to vehicle pitch rate and elastic effects, this angle is dependent upon 
the position sj along the vehicle’s lower forebody. Expressing the local velocity as a 
function of sj, an explicit expression for the local flow deflection angle 0Lf(si) is obtained. 
Similarly, the control surface deflects the air flow by an angle 0 which is dependent upon 
control surface location. (See Figure 3.) Here the control surface deflection rate 5 and its 
chord length are assumed negligible. This allows the local flow deflection angle at the 
control surface ©l cs to be considered constant for the entire surface. 


In computing the forces on the vehicle, the trigonometric sine of each flow deflection angle 
is sought rather than the angle itself. These quantities are defined as: 


sin(0 Lf <sO)=%^ sin(9 Lcs )= ^ cs ‘^ (4.2-1) 

\\lM |vj 

For the lower forebody: 

V Lf (si) = V=, + Vq(si) + Vq(si), fif = sin('ti+AxiT|)i + cos(Xi+AXiTl)k 

— ♦ A ^ 

where, V„=Vo=cos(a)i+V„sin(a)k 

V q (si)=a> rb x r(si); (Or b =qj, r(si)={ sicos(Xi+XiTO-xi }i+{ zi-sisin(xi+Xiri) }k 
Vf|(si)=-siAxiTin 

Thus the local velocity at the lower forebody VlK s 0 is: 

VLt<si)={V«cos(a)+ziq-si(q+Axifi)sin(Xi+AxiTi)}i + 
{Vo,sin(a)+xiq-si(q+AxiTi)cos(Xi+AxiTi)}k 



Figure 3. Geometric description of the hypersonic vehicle velocities, 
and sin 2 (©Lf(si)) is determined from equation (4.2-1): 
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sin^ufa)) = a2 + 2abs i +b2 s? = i+ 

c 2 +2absi+b 2 s^ c 2 +2absi+b 2 s 2 

where, a = VooSin(a-Kti+AxiT|)+q{zisin(xi+AxiT|)+xicos(xi+Axir|)} 

b = -(q+Axiq) 

c 2 = {VooCos(a)+qzi} 2 +{VooSin(a)+qxi} 2 
For the control surface: 

Vlcs = Voo + Vq + V-q, n cs =sin(S-Ax2r|)i + cos(8-Ax2Tl)k 

where, V q =co rb x r cs ; r cs =(x cs -Ax2T|r z )i+(z cs +Ax2Tlr x )k 

Vq = -Ax 2 Tlj x r^; rap=(r x -Ax 2 Tlr z )?+(r z +Ax2Tir x )k 

With these expressions the local velocity at the control surface y Lcs is computed: 


Vlcs={ VooCos(a)+q(z cs +Ax2T|r x )-Ax2Ti(r z +Ax2r|r x ) } i + 
{ V„sin(a)-q(x cs -Ax2T|r z )+Ax2Tl(r x -Ax2Tir z ) }k 

and sin 2 ( 0 L CS ) is determined from equation ( 4 . 2 - 1 ): 


sin 2 (0 Lcs ) = 


{sin(g+8-Ax2'n) + + b cs ) 2 

{c cs } 2 +{d cs } 2 


where, ac S = [(Ax2Tlr x +z cs )sin(8-Ax2Tl)+(Ax2'nr z -x cs ) cos(8-Ax2Tl)] 

V CQ 
A ’ 

b cs = TJ 1 * [zisin(8-Ax2Tl)+xicos(8-Ax2Tl)] 

V OO 

A ’ 

c cs = cos(a)+- 3 - (z cs +Ax 2 Tir x ) - — -p- (r z +Ax 2 Tir x ) 

’ OQ ’ OO 

d cs = sin(a) - (x cs -Ax 2 rir z ) + pp- (r x -Ax 2 Tlr z ) 


( 4 . 2 - 2 ) 


( 4 . 2 - 3 ) 
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To compute the aerodynamic forces and moments acting on the lower forebody, the 
pressure distribution along the lower forebody is integrated over si. In equations (4.1-1)- 
(4.1-4) two integrals appear in the aerodynamic forebody forces which reduce to: 


Hi 

sin 2 (0Lf(si))dsi 

0 


f lt 

I sisin 2 (0Lf( s i))dsi 

Jo 


These two integrals make apparent the need to express sin 2 (0Lf(si)) as an explicit function 
of sj. Equation (4.2-2) exhibits this function as a ratio of two quadratic polynomials in Sj, 
where each of the above integrals may be expressed as a sum of two separate integrals: 



Essentially, to compute the aerodynamic forces and moments on the vehicle's lower 
forebody requires analytic solutions of four integrals, two of which have trivial solutions 
while the other two are determined from an integral table. 





c 2 +2absj+b 2 sf 


si 


From [5] the closed form solutions of Ii and I 2 depend on the values of a, b and c. For 
computational purposes there are two cases that may arise. One case occurs when either 
a 2 =c 2 or the constant (c 2 ) in the integrand's denominator is much greater than the terms 
containing s\. A conservative measure for the second condition uses the largest value for 
si, i.e., c 2 » 2abli+(bli) 2 . The second case occurs when both of these conditions are 
violated: 


Case I. a 2 =c 2 or c 2 »2abli+(bli) 2 
Case II. aPtc 2 and c 2 is not » 2abli+(bli) 2 


Case I. 
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Case II. 

Ii 

i 


h 

sin 2 (0 Lf (si))ds 1 = (&f U , 



sisin 2 (0 Lf (si))dsi = ^(^-f 


let a=c 2 -a 2 . 


o>0 

j-^[tan-i(^J±^)-tan 
|b|Va |b|Va 


SZ<Q 


-u_ali 


( m )1 




b 2 li 

ab-lbjiTa 


+1 


b 2 I 


ab-^bjV^a 


J — +1 


I 2 = -l-ln{(^f+2^.1 1 + l}-Al 1 
2b 2 c c 2 b 


li 

sin 2 (0Lf<si))dsi = h - alj , 


i 


i 2 

sisin 2 (0Lf(si))dsi =-^~ - ch 


4.3. Modification of Control Surface Forces and Moments 


Expressions for the forces and moments produced by the control surface were presented in 
[2]. But these expressions do not include: l)Effect of elastic motion at the control surface 
location, 2)Contribution of the control surface force to the first elastic generalized force. 
Modifications to include these effects are derived here. 

Considering the elastic motion at the control surface location will modify the local flow 
deflection angle 0 l cs . In Section 4.2 elastic motion was considered and the resulting 

expression for sin 2 (0L CS ) includes this effect. The effect of deflection rate 8 and chord 
length of the control surface on the pressure distribution are considered negligible, and the 
local flow deflection angle is considered constant over the control surface chord. This in 
turn yields a constant pressure distribution for the entire control surface. Equations (4.1- 
5)-(4.1-8) (taken from [2]) are modified by replacing sin 2 (a+5) with sin 2 (0L CS ) of equation 
(4.2-3). Also 8, the control surface deflection, becomes (8-AT2TI) so that: 

Xcs= -q oo Cp N sin 2 (0Lcs)sin(8-Ax2Tl)^ 

Z cs = -q oo C PN sin 2 (0Lcs)cos(8-AT 2 Tl)^- 
The pitching moment produced by the control surface is: 


(4.3-1) 

(4.3-2) 


M cs - X cs (z cs +Al2tir x ) ■ Z cs (x ts -Ax 2 Tir z ) (4.3-3) 

The control surface contribution to the elastic generalized force was not considered in [2]. 
The virtual work produced by the control surface force can be differentiated with respect to 
the generalized displacement to yield the generalized force Qc S contributed by the control 
surface: 

^ _38W d(F cs -Ar) 

dr | dr| 

— ► A /V 

where, F^s— -t* Zc§k 

Ar = -AT 2 T|r z i + Ax 2 "nr x k 


Qcs= - Ax 2 T|r z - Ax 2 r z X cs + -^1 Ax 2 Tir x + Ax 2 r z Z cs 


so that: 
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Differentiating X cs and with respect to ti yields: 

= q~Cpn^ [Ax 2 co S (5-Ax 2 T)) S in 2 (e Lcs ) - sin(5-Ax 2 T O 3sin ^ 0Lcs) ] 

= - qooCpn^- [Ax 2 sin(6-Ax 2 ri)sin 2 (0 Lcs ) - cos(5-Ax 2 il) 3sin ^ 9Lcs ^ ] 

Differentiating equation (4.2-3) with respect to 11 and substituting back into the above 
equations results in a lengthy expression for Q cs : 

(Ax 2 r|r x -r z )X cs + (Ax 2 T)r z +r x )Z cs 

, oa- ^ _ v x r {q(xi+Ax 2 rir z )-Ax2TlAx 2 f|r z }sin(5-Ax 2 Ti) 

+ zax2T\ (r x z, cs -r z A cs [ 

V^(sin (a+S-Ax^l+acs+bcs) 

| {Ax 2 qAx 2 T)r x -q(zi+Ax 2 Tir x )}cos(5-Ax 2 Ti) 

V 00 (sin(a+5-Ax 2 ii)+a cs +b cs ) 

V 00 cos(a+5-Ax 2 q) 

V^CsinCa+S-Ax^J+acs+bcs) 

(q- A x 2 T) )(r x c cs +r z d cs ) 

V«.(Ccs + dcs) 





5. COMPUTING THE AEROPROPULSIVE FORCES AND MOMENTS 


5.1. Computation of Useful Variables 

The aeropropulsive force equations (4.1-1) to (4.1-16) require various quantities which 
depend on the vehicle states, controls and physical parameters. Many of these quantities 
are computed in the block COMPVAR. For example, the freestream temperature and 
pressure are functions of altitude and are determined according to the ARDC Standard 
Atmosphere. Also determined are freestream velocity, Mach number and dynamic 
pressure, as well as angle of attack. In addition, COMPVAR calculates the instantaneous 
elastic deflection of the vehicle and its elastic deflection rate from the first generalized elastic 
co-ordinate and its mode shape. 

5.2. Aerodynamic Body Forces and Moments 

As noted in Chapter 4, the aerodynamic forces acting on a hypersonic vehicle are expressed 
in terms of the local flow deflection angle 0 along the forebody and at the control surface. 
The block THETAL determines sin 2 (0Lf) at the forebody by computing the coefficients a, b 
and c in equation (4.2-2). It then computes the closed-form solutions for the integral 
terms, as discussed in Subsection 4.2.2. The block AEROB utilizes equations (4.1-1) to 
(4.1-4) to compute the aerodynamic forces and moments produced by the pressure 
distribution on the vehicle's lower forebody and upper surface. 

5.3 Aerodynamic Control Surface Forces and Moments 

The block AEROC calculates sin 2 (0 Lcs ), using the parameters a^, b cs , c cs , and d cs in 
equation (4.2-3), and computes the control surface forces and moments in accordance with 
equations (4.3-1) to (4.3-4). 


5.4 Engine Thermodynamics and Thrust Forces and Moments 

The block ENGTHERM is made up of four sub-blocks: STAGE1, STAGE2, STAGE3, 
and ST AGEE. Given freestream conditions, engine control settings (A D and T 0 ), and the 
local flow deflection angle at the engine inlet (0Lf<si=O)), the functions in ENGTHERM are 
used to determine inlet conditions at the engine diffuser, combustor, nozzle, and at the 
engine exit. The equations (5.4-1) to (5.1-12) below are used [2]. A numerical search 



using Newton’s method is used to solve equations (5.4-4) and (5.4-10) iteratively for the 
Mach number. The block THRUST then computes the engine-module thrust force and 
moment using equations (4.1-9) to (4.1-12). 

External Diffuser: 


Mi = — 

M oo cos(0 Lf ) 

: (5.4-1) 

V 

1+2 (Y'l)Mo. 2 sin 2 (0 L f) 

Pi = p a 

| 1+ YCpMoo 2 sin 2 (0Lf) 

(5.4-2) 

Tl = Poo 

’ l+l(Y-l)M3in 2 (0 L f) 

(5.4-3) 


Internal Diffuser: 


Combustor: 


y+ i Y+ i 

(l+i(Y-l)M 2 2 )^T__ 2 (l+^(Y-l)M! 2 )^r 


M 2 2 
P 2 = Pi 


1+lfY-DM! 2 


M] 2 

tfr 


i+1(y-dm 2 2 


T 2 = Tj 


l+ifY-DMi 2 


l+i(Y-l)M 2 2 


M 3 2 (l4-(Y-1)M 3 2 ) _ M^l+lfY-DM^) _M 2 2_ Ia 
( 1 +YM 3 2 f ( 1 +yM 2 2 f ( 1 +yM 2 2 ) 2 T 2 


P 3 =Pd 


i+ym 2 2 

i+ym 3 2 _ 


(5.4-4) 


(5.4-5) 


(5.4-6) 


(5.4-7) 


(5.4-8) 


^ ^[(i+ym 2 2 ) m 3 ? 

T3 T {(l^i7)Mlj 


(5.4-9) 



Internal Nozzle: 



Pe = P3 


l+^-(Y-l)M 3 2 

1+1(Y-1)M 2 2 


_Y_ 

Y-l 


T e = Td 


1+1(Y-1)M 3 2 

l+lCY-DMe 2 


(5.4-10) 

(5.4-11) 


(5.4-12) 


5.5. External Nozzle Forces and Moments 

The engine exhaust plume is bounded above by the vehicle lower aftbody, and below by a 
free shear layer. This variable structure is termed the external nozzle. The pressure 
distribution along the lower aftbody produces forces and moments on the vehicle aftbody. 
The block NOZZLE computes these quantities using equations (4.1-13) to (4.1-16), taken 
from [2]. 


6. VEHICLE CONFIGURATION AND FLIGHT CONDITION 


6.1. Model Parameters 

This section lists several pertinent constants associated with the spherical rotating earth and 
ARDC standard atmosphere models. Many atmospheric parameters are dependent on 
altitude and are contained in the function strdatm. Also given are vehicle dimensions, 
elastic and mass properties, and engine parameters for a particular vehicle configuration. 
All of the values below are defined in the script file inpar. 


Spherical Rotating Earth: 


Symbol 
Text Code 

Value 

Description 

Units 

Re 

re 

2.09256E+6 

Earth radius 

[ft] 

He 

mu 

1.40764E+16 

Earth gravitational parameter 

[ft 3 /s 2 ] 

Ofc 

ome 

7.297205E-5 

Earth rotation rate about its axis 

[rad/s] 


ARDC Standard Atmosphere (1959): 


Symbol 
Text Code 

Value 

Description 

Units 

Y 

gam 

1.43 

Ratio of specific heats 

[-] 

Ra 

gasc 

1.716545E+3 

Atmospheric gas constant 

[ft 2 /°R-s 2 ] 

Cpn 

cpn 

2.0 

Pressure coefficient 

[-] 


(A complete model for a spherical rotating earth and ARDC Standard Atmosphere is given 
in [1].) 



Vehicle Dimensions and Mass Properties (Example): 


Symbol 
Text Code 

Value 

Description 

Units 

*1 

taul 

14.0 

Vehicle nose angle 

[deg] 

*2 

tau2 

20.0 

Vehicle tail angle 

[deg] 

Axi deltaul 

1.0 

Forebody elastic mode shape 

[deg] 

AT2 

deltau2 

1.0 

Aftbody elastic mode shape 

[deg] 

L 

1 

150 

Total vehicle length 

[ft] 

h 

h 

22.2 

Total vehicle height 

[ft] 

Li 

11 

89.03 

Vehicle forebody length 

[ft] 

l 2 

12 

60.97 

Vehicle aftbody length 

[ft] 

Xcg 

xcg 

90.0 

c.g. position aft of vehicle nose 

[ft] 

Zcg 

zcg 

11.25 

c.g. position below vehicle nose 

[ft] 

xcs 

xcs 

-52.5 

Control surface position w.r.t. center 
of gravity in Body x-axis 

[ft] 

Zcs 

zcs 

-11.25 

Control surface position w.r.t. center 
of gravity in Body z-axis 

[ft] 

Scs 

scs 

22.5 

Control surface area (per vehicle width) [ft 2 /ft] 

m 

m 

500.0 

Vehicle mass (per vehicle width) 

[slug/ft] 

lyy 

iyy 

1.0E+6 

Vehicle Body y-axis moment of inertia 

(per vehicle width) [ft 2 -slug/ft] 

m£i 

mel 

40.0 

Generalized mass associated with the 
1 st generalized elastic DoF 
(per vehicle width) 

[slug/ft] 

coi 

oml 

18.0 

Undamped natural frequency of 
1 st generalized elastic DoF 

[rad/s 2 ] 

Cl 

zetl 

0.01 

Damping ratio of 1 st generalized 
elastic DoF 

[-] 

Engine Parameters (Example): 



Symbol 
Text Code 

Value 

Description 

Units 

An 

an 

6.35 

Nozzle area ratio 

[-] 

Ae 

ae 

8.88 

Nozzle exit area (per vehicle width) 

[ft 2 /ft] 


(A complete description of the geometry for this particular vehicle is given in [2].) 


6.2. Initial Conditions-Example 
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The formulation of ASUHS1 is based upon a wings level (<|>=0°) flight condition. For the 
example presented here, the initial vehicle flight-path angle is zero ( y=0 ° , 0=a), with the 
velocity- vector directed eastward (\|/=90°). This information specifies the vehicle's initial 
orientation. 

Furthermore, suppose the vehicle position is specified with an altitude of 85,000 ft., over 
the intersection of the equator (latitude X=0°) and prime meridian (longitude x=0°). (See 
Figure 2.) This information specifies the three vehicle states h, X and x. The vehicle's 
initial Mach number for this flight condition is 8.0. The ambient conditions according to 
the ARDC Standard Atmosphere at are: 

P M = 45.82 lb/ft 2 
^ = 394.3 °R 

Moo = 8.0 

Voo = 7870.5 ft/s 

At this point a trimmed flight condition (operating point) must be determined. Appendix C 
presents a trimming method which may be used to determine this trimmed flight condition. 
The method assumes a particular flight condition structure. Here the flight condition 
structure is defined as the statement of those variables that are specified ("fixed") and those 
that are allowed to vary ("free"). The flight condition structure used for this method 
specifies: vehicle position (h, x, A.), Mach number (M^,), roll and yaw attitude (<J>, \}/), 
flight path angle (y), and engine temperature control (T 0 ) and nozzle area ratio (An). The 
vehicle rotation and elastic rates (q, T|), and accelerations (li, w, q, rj) are specified to be 
zero. The chosen free variables are: angle of attack (a), pitch control surface deflection 
(8), diffuser area control (Ad), and elastic deflection (T|). Under this scenario the resulting 
trimmed vehicle states and controls are: 

h G = 85,000 [ft] 
u 0 = 7,806 [ft/s] 

w 0 = -1002 [ft/s] (a 0 =-7.317°) 
q 0 = 0.0 [rad/s] 
llo= 1.243 [-] 

f | 0 = 0.0 

Xq = 0.0 [rad] 
x 0 = 0.0 [rad] 


Pl 0 = 0.04512 [-] 

P 2o = -0.04512 [-] 

P 3o = 0.7057 [-] 
p 4o = 0.7057 [-] 

S 0 = 25.21 [deg] 

A Do = 0.5004 [-] 

T 0 = 2000 [R°] 

Due to round-off of the free variables, the accelerations do not exactly result in zero values. 
The following accelerations result from the given initial conditions and controls: 

Uo = 4.475E-4 [ft/s 2 ] 
w 0 = 7.514E-4 [ft/s 2 ] 
q 0 = 1.336E-5 [rad/s 2 ] 
ii 0 = -1.402E-2 [1/s 2 ] 

The script file initstate. is currently set up to perform an example simulation using the 
above initial conditions and a doublet input for 5 (shown below in Figure 4). The resulting 

trajectory is presented in Figure 5. 
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Figure 4. Doublet input for Control Surface Deflection Angle 8 (±1 degree from trim). 
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Figure 5. Simulation time responses to a doublet input in 8. 








7. DISCUSSION AND POTENTIAL PROBLEMS 
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The expressions derived in this work rely entirely upon the aeroelastic model presented in 
Reference [2]. This elastic model was intended to provide the analyst with first-order 
effects of the expected interactions between the engine and aircraft and also between the 
rigid body and elastic degrees of freedom. Thus, if the elastic model is changed greatly, 
the present model may or may not capture all first-order effects. Certainly, if the elastic 
model is altered, the expressions derived herein may no longer be valid. 

A major contribution of this work has been to integrate the pressure distribution along the 
vehicle's lower forebody. It was shown that the integrals of the pressure distribution 
produced by this hypersonic vehicle model have closed-form solutions. The structure of 
these solutions depends upon the vehicle's dynamic state (Case I or Case II in Subsection 
4.2.2). In computing these solutions, the magnitudes of some parameters may be small but 
non-zero. A measure for zero (to be adjusted at the user's discretion) must be 
incorporated. From experience, a term smaller than IE-6 (the current zero measure in the 
function intheta.) will be considered as zero. The current zero measure provides a smooth 
numerical transition from Case I to Case II. 

The derivation of the engine thermodynamic conditions [2] relies on the a priori assumption 
of non-choked flow (i.e., Mach numbers always greater than unity throughout the engine- 
module). This certainly is the case for the vehicle configuration presented here, but in 
general non-choked flow is not guaranteed for arbitrary engine control settings. This issue 
may become critical when attempting to trim about a new flight condition. In addition, if 
the control settings are allowed to vary (e.g., in a closed-loop control simulation), the 
analyst must ensure that engine controls do not choke the fluid flow and thus violate the 
engine modeling assumptions. If the engine controls do choke the fluid flow, a "choked 
flow" warning is printed to MATLAB®'s command window. If this occurs it is suggested 
that the simulation be cancelled by typing "control-C," if the simulation does not cease by 
itself. 

One peculiarity encountered while running ASUHS1 is the initial warning that is flagged to 
MATLAB®'s command window. This warning is a result of using the previous 
integration time step values of M 2 and M e as the initial guesses in the Newton iteration 
search. This method reduces computation time but at the same time introduces this 
annoying warning. This warning, due to this implementation of the Newton iteration 
search, may be ignored. Due to this implementation in stage2 and stagee, the simulation 
start time must be set at zero. 
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APPENDIX 


A. DESCRIPTION OF ASUHS1 MATLAB® FUNCTIONS 
Function Name: strdatm 
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Description: This function implements the 1959 ARDC Standard Atmosphere to 

compute the freestream temperature and pressure given the vehicle’s 
altitude. This model was obtained from Reference [1]. 


Input Data: 

Symbol Description Units 

Text Code 

h h Geometric altitude as measured from the earth’s surface [ft] 

R e re Earth’s radius [ft] 

g 0 go Gravitational acceleration at the earth’s surface [ft/s 2 ] 

Ra gasc Atmospheric gas constant [ft 2 /°R-s 2 ] 


Output Data: 

Symbol Description 

Text Code 

Too tinf Atmospheric freestream temperature 

Poo pinf Atmospheric freestream pressure 


Units 

[°R] 

[lb/ft 2 ] 


Procedure: Geopotential altitude hg is first computed from geometric altitude: 


h ‘ = h <R^> 


With a geopotential altitude the atmospheric layer is determined and the 
corresponding values for T re f, C, h re f and P re f are determined according 
the ARDC Standard Atmosphere. These four parameters are used to 
compute the freestream temperature and pressure: 


Too — T re f + C(hg-h re f) 


for C=0, 
for C*0, 


Poo=Pref {exp[- (hg-href)] } 

Poo=P re fj 1 + (hg-h re f) (crJ 


Function Name: speedalpha 
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Description: This function computes the total vehicle velocity and angle of attack 

given the vehicle velocity components. 

Input Data: 

Symbol Description Units 

Text Code 

u u Vehicle forward speed [ft/s] 

w w Vehicle sink rate [ft/s] 

Output Data: 

Symbol Description Units 

Text Code 

Voo vinf Total vehicle speed (in pitch plane) [ft/s] 

a al Vehicle angle of attack [rad] 

Procedure: The total vehicle speed is computed as: 

Voo=Vu 2 +w 2 

The vehicle angle of attack is determined by: 

a=tan’ 1 (^) 



Function Name: machq 
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Description: This function computes the freestream Mach number and freestream 

dynamic pressure. 


Input Data: 

Symbol 

Description 

Units 

Text Code 

Tm tinf 

Atmospheric freestream temperature 

[°R] 

Poo pinf 

Atmospheric freestream pressure 

[Ib/ft 2 ] 

Voo vinf 

Total vehicle speed (in pitch plane) 

[ft/s] 

Ra gasc 

Atmospheric gas constant 

[ft 2 /°R-s 2 ] 

Y gam 

Ratio of specific heats 

[-] 

Output Data: 

Symbol 

Description 

Units 

Text Code 

Moo minf 

Freestream Mach number 

[-] 

qoo qinf 

Freestream dynamic pressure 

[Ib/ft 2 ] 


Procedure: The speed of sound asonic is computed using the isentropic flow equation: 

a sonic = ^ YRaT^ 

With the sonic speed the freestream Mach number is determined from: 

M V “ 

Moo= 

a sonic 

Freestream dynamic pressure is then computed from: 

Qm = 2 P. coM^ 



Function Name: elasdef 


Description: This function computes the elastic deflections and deflection rates for 

the forebody and aftbody sections. 


Input Data: 


Symbol 
Text Cods 

Description 

Units 

*n 

eta 

First generalized elastic co-ordinate 

[-] 


deta 

First generalized elastic co-ordinate rate 

[1/s] 

Axi 

deltaul 

Forebody elastic mode shape 

[rad] 

AX2 

deltau2 

Aftbody elastic mode shape 

[rad] 

Xi 

taul 

Vehicle nose angle 

[rad] 

*2 

tau2 

Vehicle tail angle 

[rad] 

Output Data: 



Symbol 
Text Code 

Description 

Units 

AXlT] 

deltaunl 

Elastic deflection of forebody 

[rad] 

AX 2 TI 

deltaun2 

Elastic deflection of aftbody 

[rad] 

'tn 

ttaul 

Total (effective) vehicle nose angle 

[rad] 

Xt 

ttau2 

Total (effective) vehicle tail angle 

[rad] 

Axii) 

dtaul 

Elastic deflection rate of forebody 

[rad/s] 

AX2T1 

dtau2 

Elastic deflection rate of aftbody 

[rad/s] 


Procedure: The elastic deflections are computed by multiplying the mode shape 

components with the elastic co-ordinate, where the first component 
corresponds to the forebody section and second component corresponds 
to the aftbody section. The elastic deflections are then added to the 
nose and tail angles to yield the "effective" nose and tail angles of the 
vehicle: 


X n = Ti + ATlT| 

X t = X2 + AX2T1 

The elastic deflection rates are computed by multiplying the 
corresponding mode shape components with the elastic co-ordinate rate. 


Function Name: coef 
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Description: This function computes the "coefficients" involved in the expression for 

sin 2 (9 L (si)). 


Input Data: 

Symbol Description Units 

Text Code 

V«. vinf Freestream velocity [ft/s] 

a al Vehicle angle of attack [rad] 

T n ttaul Total (effective) vehicle nose angle [rad] 

Axiii dtaul Elastic deflection rate of forebody [rad/s] 

q q Vehicle pitch rate [rad/s] 

xi xl Location of c.g. with respect to lower apex in body x-axis [ft] 

z\ zl Location of c.g. with respect to lower apex in body z-axis [ft] 

Output Data: 

Symbol Description Units 

Text Code 

a a Represents the local flow normal to lower forebody at the [-] 

engine inlet 

b b Represents the effective pitch rate along the lower forebody [1/ft] 

surface 

c 2 csqrd Represents the total velocity squared at the engine inlet [ 1/ft 2 ] 


Procedure: The three "coefficients" are computed according to equation (4.2-2) in 

Section 4.2.1. 



Function Name: intheta 
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Description: This function computes the algebraic solutions of the integrals involving 

sin 2 (0 L (si)). 


Input Data: 

Symbol Description Units 

Text Code 

a a Represents the local flow normal to lower forebody at the [-] 

engine inlet 

b b Represents the effective pitch rate along the lower forebody [1/ft] 

surface 

c 2 csqrd Represents the total velocity squared at the engine inlet [ 1/ft 2 ] 

ll ell Length along vehicle forebody [ft] 

Output Data: 

Symbol Description Units 

Text Code 

sqth sin 2 (0 L (si=O)) i.e., computed at the engine inlet [-] 

isqth Integral of sin 2 (0 L (s]))dsj from sj=0 to S|=l^ [ft] 

issqth Integral of Sisin 2 (0L(si))dsj from sj=0 to s |=lj [ft 2 ] 


Procedure: This function begins by first setting an acceptable measure for zero, (Was 

(currently 0 me as = lE-6). The parameters a and s are then computed as: 


a = c 2 - a 2 
s = V|oj 


The logic to determine which algebraic solution applies is as follows: 


If s > Omeas and c 2 < 


2abli+(bli)' 


0 , 


meas 


If C7> 0 
Case II (o>0) 

Else 

Case II (a<0) 

Else 


Case I 



Comparing s (rather than a) to 0 mea s permits use of the same 0 meas value 
to be used in the two inequalities in the first line of this logic sequence. 
The implementation of the analytic solutions in Section 4.2.2 avoids, 
whenever possible, division by Ibl, a or Va since these quantities are 
typically very small. 

After computing the appropriate algebraic solutions the term sin 2 (0 Lf <O)) 
is computed as: 

sin 2 (0Lf(si= 0 )) = ai 
c 2 
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Function Name: stagel 

Description: This function computes the thermodynamic conditions (Mi, Pi, Ti) at the 

engine inlet (diffuser inlet) referred to as Stage 1. Newtonian Impact 
Theory was applied in the derivation of these conditions [2]. 


Input Data: 


Symbol 
Text Code 

Description 

Units 

qoo 

qinf 

Freestream dynamic pressure 

[lb/ft 2 ] 

Mco 

minf 

Freestream Mach number 

[-] 

Poo 

pinf 

Atmospheric freestream pressure 

[lb/ft 2 ] 

Too 

tinf 

Atmospheric freestream temperature 

[•R] 

- 

sqth 

Evaluation of sin^Gjjsj)) at sj=0 (engine inlet) 

[-] 

Y 

gam 

Ratio of specific heats 

[-1 

Cpn 

cpn 

Pressure coefficient 

N 

Output Data: 



Symbol 
Text Code 

Description 

Units 

Mi 

ml 

Mach number at engine inlet (diffuser inlet) 

[-] 

Pi 

pl 

Pressure at engine inlet (diffuser inlet) 

[lb/ft 2 ] 

Ti 

tl 

Temperature at engine inlet (diffuser inlet) 

[*R] 


Procedure: The function first determines if the external diffuser has choked the fluid 

flow at the engine inlet. If choked flow does occur then a warning is 
flagged to the command window. 

To facilitate computation a useful parameter Mgoo is computed: 

Mgoo = 1 + - ML 

Equations (5.4-l)-(5.4-3) in Section 5.4 are then implemented to compute 
Mi, Pi and Ti, respectively. 


Function Name: stage2 
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Description: This function computes the thermodynamic conditions (M 2 , P 2 , T 2 ) at the 

exit to the internal diffuser (combustor inlet) referred to as Stage 2. A 
supersonic isentropic flow assumption was made in the derivation of 
these conditions [2]. To solve for M 2 a Newton iteration search was 
implemented. 

Input Data: 

Symbol Description Units 

m . y— 1 1 


Text 

Ad 

Code 

ad 

Diffuser area ratio control 

[-] 

Mi 

ml 

Mach number at engine inlet (diffuser inlet) 

[-] 

Pi 

Pi 

Pressure at engine inlet (diffuser inlet) 

[lb/ft 2 ] 

Ti 

tl 

Temperature at engine inlet (diffuser inlet) 

[°R] 

Y 

gam 

Ratio of specific heats 

[-] 

Ym 

gamm 

A useful constant for computation, y m ==y/(Y- 1 ) 

[-] 

Ys 

sph 

A useful constant for computation, y s =(y-i- 1 )/(y- 1 ) 

[-] 

t 

time 

Simulation time 

[sec] 

M 2i 

m2i 

First initial guess for the Newton iteration search 

[-] 

Output Data: 

Symbol 

Description 

Units 

Text 

Ad 

Code 

ad 

Diffuser area ratio control 

[-] 

m 2 

m2 

Mach number at engine stage 2 (diffuser exit) 

[-] 

P2 

p2 

Pressure at engine stage 2 (diffuser exit) 

[lb/ft 2 ] 

t 2 

t2 

Temperature at engine stage 2 (diffuser exit) 

m 


Procedure: The function first computes the useful parameter Mg 1 : 

Mgi = 1+ ~-Mi 

It then determines if the diffuser area ratio control will choke the fluid 
flow at the diffuser exit. If choked flow occurs a warning is flagged to 
MATLAB®'s command window. If the value for Ad maintains 
supersonic flow, a Newton iteration search is utilized to solve for M 2 . 



To reduce the number of cycles executed by the Newton iterative search 
the previous value for M 2 is used as the initial guess to start the search. 
To overcome the problem during the initial time calculations (there is no 
previous M 2 available), the following logic sequence is used: 


If t >0.0 


M 2 = M2i 


Else 


M 2 = 


Mi 

2.5 


The function then goes on to iteratively solve for M 2 using equation (5.4- 
4) in Section 5.4, via Newton's method. 


With a value for M 2 the useful parameter Mg 2 is computed: 

Y-l 2 

Mg 2 = 

Equations (5.4-5) and (5.4-6) in Section 5.4 are then implemented which 
compute P 2 and T 2 , respectively. 


Function Name: stage3 
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Description: This function computes the thermodynamic conditions (M 3 , P 3 , T 3 ) at the 

exit to the combustor (nozzle inlet) referred to as Stage 3. A constant 
area heat addition assumption was made in the derivation of these 
conditions [2], where T 0 is the total temperature change across the 
combustor. M 3 is determined via the quadratic formula. 


Input Data: 

Symbol 

Description 

Units 

Text 

T 0 

Code 

to 

Total temperature change control (across combustor) 

[°R] 

m 2 

m 2 

Mach number at engine stage 2 (combustor inlet) 

[-1 

P2 

p 2 

Pressure at engine stage 2 (combustor inlet) 

[lb/ft 2 ] 

t 2 

t 2 

Temperature at engine stage 2 (combustor inlet) 

[°R] 

y 

gam 

Ratio of specific heats 

[-] 

Output Data: 

Symbol 

Description 

Units 

Text 

T 0 

Code 

to 

Total temperature change control (across combustor) 

ra 

m 3 

m3 

Mach number at engine stage 3 (combustor exit) 

[-] 

P 3 

p3 

Pressure at engine stage 3 (combustor exit) 

[lb/ft 2 ] 

t 3 

t3 

Temperature at engine stage 3 (combustor exit) 

[°R] 


Procedure: The function begins by computing useful parameters, g2 and Mg 2 , which 

are used continuously throughout this routine: 

g 2 = 1+ yM 2 , Mg2 _ 1+ 2 

The total temperature change control is checked to verify that the flow 
within the combustor is non-choked. If choked flow occurs, a warning is 
flagged to MATLAB®'s command window. If the value for T 0 maintains 
supersonic flow, then equation (5.4-7) in Section 5.4 is represented as a 
quadratic function in M 3 and the quadratic formula applied to solve for 
the appropriate root, i.e., the real positive root which is greater than unity. 



With a value for M3 the useful parameter g3 is computed: 
g 3 = I+7M3 


The equations (5.4-8) & (5.4-9) in Section 5.4 are then implemented 
compute P3 and T3, respectively. 


Function Name: stagee 

Description: This function computes the thermodynamic conditions (M e , Pe, T e ) at the 

exit to the internal nozzle (engine exit) referred to as Stage e. A 
supersonic isentropic flow assumption was made in the derivation of 
these conditions [2], To solve for M e a Newton iteration search was 
implemented. 


Input Data: 

Symbol Description Units 

Text Code 


An 

an 

Nozzle area ratio 

[-] 

m 3 

m3 

Mach number at engine stage 3 (nozzle inlet) 

[-] 

P3 

p3 

Pressure at engine stage 3 (nozzle inlet) 

[lb/ft 2 ] 

t 3 

t3 

Temperature at engine stage 3 (nozzle inlet) 

[°R] 

y 

gam 

Ratio of specific heats 

[-1 

Ym 

gamm 

A useful constant for computation, y m =y/(Y-l) 

[-] 

Ys 

sph 

A useful constant for computation, y s =(Y+l)/(y-l) 

[-] 

t 

time 

Simulation time 

[sec] 

M ei 

mei 

First initial guess for the Newton iteration search 

[-] 

Output Data: 

Symbol 

Description 

Units 

Text 

An 

Code 

an 

Nozzle area ratio 

[-] 

M e 

mex 

Mach number at engine exit (nozzle exit) 

[-] 

Pe 

pex 

Pressure at engine exit (nozzle exit) 

[lb/ft 2 ] 

T e 

tex 

Temperature at engine exit (nozzle exit) 

[°R] 


Procedure: The function first computes the useful parameter Mg 3 : 

Mg 3 =1+^M^ 

It then determines if the nozzle area ratio will choke the fluid flow at the 
nozzle exit. If choked flow occurs a warning is flagged to MATLAB®'s 
command window. If the value for Ajsj maintains supersonic flow a 
Newton iteration search is utilized to solve for M e . 



To reduce the number of cycles executed by the Newton iterative search 
the previous value for M e is used as the initial guess to start the search. 
To overcome the problem during the initial time calculations (there is no 
previous M e available) the integration time variable t is also used as input 
to this function. The logic sequence is as follows: 

If t >0.0 

Mg — Mei 

Else 

Me = 3M 3 

The function then goes on to iteratively solve for M e using equation (5.4- 
10) in Section 5.4, via Newton's method. 

With a value for M e the useful parameter Mg e is computed: 

Mge =1+^M| 

The equations (5.4-1 1) & (5.4-12) in Section 5.4 are then implemented to 
compute P e and T e , respectively. 


Function Name: aerobody 
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Description: This function computes the aerodynamic forces and moments acting on 

the vehicle forebody. Note these forces are determined per unit width of 
the vehicle. 


Input Data: 



Symbol 
Text Code 

Description 

Units 

Qoo 

qinf 

Freestream dynamic pressure 

[lb/ft 2 ] 

P=o 

pinf 

Atmospheric freestream pressure 

[lb/ft 2 ] 


ttaul 

Total (effective) vehicle nose angle 

[rad] 

- 

isqth 

Integral of sin 2 (0L(si))dsj from si=0 to Sj=l| 

[ft] 

- 

issqth 

Integral of Sisin 2 (0L(si))dsj from sj=0 to Sj=I| 

[ft 2 ] 

ll 

ell 

Length along vehicle forebody 

[ft] 

Cpn 

cpn 

Pressure coefficient 

[-] 

xi 

xl 

Location of c.g. with respect to lower apex in body x-axis 

[ft] 

zi 

zl 

Location of c.g. with respect to lower apex in body z-axis 

[ft] 

L 

1 

Vehicle length 

[ft] 

Xcg 

xcg 

Distance from vehicle nose to c.g. 

[ft] 

Axi deltaul 

Forebody elastic mode shape 

[rad] 

Output Data: 



Symbol 
Text Code 

Description 

Units 

X A 

xa 

Body x-axis aerodynamic force acting on vehicle forebody 
(per unit width of the vehicle) 

[lb/ft] 

Za 

za 

Body z-axis aerodynamic force acting on vehicle forebody 
(per unit width of the vehicle) 

[lb/ft] 

m a 

ma 

Body y-axis aerodynamic moment acting on vehicle forebody [ft- lb/ft] 
(per unit width of the vehicle) 

Qa 

qa 

First generalized force due to aerodynamic force acting on 
vehicle forebody (per unit width of the vehicle) 

[ft- lb/ft] 


Procedure: The function computes the magnitude of the total aerodynamic force 

acting on the vehicle forebody: 

( h 

F a = P 0 J 1 + qooCpn 1 sin 2 (0Lf(si))dsi 

Jo 
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The function then goes on to compute Xa, Za, Ma and Qa using 
equations (4.1-1) to (4.1-4) in Section 4.1. 


Function Name: aerocont 
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Description: This function computes the aerodynamic forces and moments acting on 

the control surface. Note these forces are determined per unit width of 
the vehicle. 


Input Data: 



Symbol 
Text Code 

Description 

Units 

qoo 

qinf 

Freestream dynamic pressure 

[Ib/ft 2 ] 

Vcc 

vinf 

Freestream velocity 

[ft/s] 

q 

q 

Vehicle pitch rate 

[rad/s] 

Ax 2 i) 

dtau2 

Elastic deflection rate of aftbody 

[rad/s] 

Ax 2 t] deltaun2 

Elastic deflection of aftbody 

[rad] 

a 

al 

Vehicle angle of attack 

[rad] 

8 

delta 

Pitch control surface deflection 

[rad] 

II 

ell 

Length along vehicle forebody 

[ft] 

Cpn 

cpn 

Pressure coefficient 

[-] 

r x 

rx 

Location of control surface w.r.t. lower apex in body x-axis 

[ft] 

r z 

rz 

Location of control surface w.r.t. lower apex in body z-axis 

[ft] 

h 

h 

Vehicle height 

[ft] 

x cg 

xcg 

Distance from vehicle nose to c.g. in body x-axis 

[ft] 

Zcg 

ZCg 

Distance from vehicle nose to c.g. in body z-axis 

[ft] 

Xcs 

xcs 

Location of control surface with respect to c.g. in body x-axis [ft] 

Zcs 

zes 

Location of control surface with respect to c.g. in body x-axis [ft] 

S CS 

SCS 

Control surface area (per unit width of the vehicle) 

[ft 2 /ft] 

Ax 2 

deltau2 

Aftbody elastic mode shape 

[rad] 

Output Data: 



Symbol 
Text Code 

Description 

Units 

Xc 

xc 

Body x-axis aerodynamic force acting on control surface 
(per unit width of the vehicle) 

[Ib/ft] 

Zc 

zc 

Body z-axis aerodynamic force acting on control surface 
(per unit width of the vehicle) 

[Ib/ft] 

M C 

me 

Body y-axis aerodynamic moment acting on control surface 
(per unit width of the vehicle) 

[ft- Ib/ft] 

Qc 

qc 

First generalized force due to aerodynamic force acting 
on control surface (per unit width of the vehicle) 

[ft- lb/ft] 


Procedure: The function calculates the individual terms of equation (4.2-2) in 

Section 4.2 in order to compute sin 2 (0L CS ). Equations (4.3-1) to (4.3-4) 
are then implemented to determine Xc, Zc, Me and Qc, respectively. 


Function Name: engthrust 
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Description: This function computes the forces and moments generated by the engine 

module thrust. It is assumed that the resultant thrust force is in the body 
x-axis direction. Note these forces are determined per unit width of the 
vehicle. 


Input Data: 


Symbol 
Text Code 

Description 

Units 

Ad 

ad 

Diffuser area ratio control 

[-] 

q«> 

qinf 

Freestream dynamic pressure 

[lb/ft 2 ] 

Mi 

ml 

Mach number at engine inlet 

[-] 

Pi 

Pi 

Pressure at engine inlet 

[lb/ft 2 ] 

M e 

mex 

Mach number at engine exit 

[-] 

Pe 

pex 

Pressure at engine exit 

[lb/ft 2 ] 

A e 

ae 

Nozzle exit area 

[ft 2 /ft] 

Y 

gam 

Ratio of specific heats 

[-] 

An 

an 

Nozzle area ratio 

[-] 

zi 

zl 

Location of c.g. with respect to lower apex in body z-axis 

[ft] 

Output Data: 



Symbol 
Text Code 

Description 

Units 

Xt 

xt 

Body x-axis engine thrust force (per unit width of the vehicle) 

[lb/ft] 

Zt 

zt 

Body z-axis engine thrust force (per unit width of the vehicle) 

[lb/ft] 

Mt 

mt 

Body y-axis engine thrust pitching moment (per unit width of [ft- lb/ft] 
the vehicle) 

Qt 

qt 

First generalized force due to engine thrust (per unit width of 
the vehicle) 

[ft- lb/ft] 


Procedure: The function computes the engine thrust forces and moments by direct 

implementation of equations (4.1-9) to (4.1-12) in Section 4.1. 



Function Name: ex tnozzle 


Description: This function computes the aeropropulsive forces and moments on the 

lower aftbody of the vehicle considered to be the upper surface of the 
external nozzle. A pressure distribution which varies as 1/(1 +S 2 ) has been 
assumed [2], where S2 is the position along the lower aftbody from the 
lower apex. Note these forces are determined per unit width of the 
vehicle. 


Input Data: 

Symbol 
Text Code 

Tt ttau2 

Poo pinf 

xi xl 

z\ zl 

P e pex 

12 el2 

%2 tau2 

At 2 deltau2 

At 2 T| deltaun2 


Description 

Total (effective) vehicle tail angle 
Atmospheric freestream pressure 

Location of c.g. with respect to lower apex in body x-axis 
Location of c.g. with respect to lower apex in body z-axis 
Pressure at engine exit (nozzle exit) 

Length along vehicle aftbody 
Vehicle tail angle 
Aftbody elastic mode shape 
Elastic deflection of aftbody 


Units 

[rad] 

[lb/ft 2 ] 

[ft] 

[ft] 

[lb/ft 2 ] 

[ft] 

[rad] 

[rad] 

[rad] 


Output Data: 

Symbol 
Text Code 

Xe xe 
Ze ze 
Me me 
Qe qe 


Description Units 

Body x-axis aeropropulsive force acting on the external [lb/ft] 
nozzle (per unit width of the vehicle) 

Body z-axis aeropropulsive force acting on the external [lb/ft] 
nozzle (per unit width of the vehicle) 

Body y-axis aeropropulsive pitching moment acting on the [ft- lb/ft] 
external nozzle (per unit width of the vehicle) 

First generalized aeropropulsive force acting on the [ft- lb/ft] 

external nozzle (per unit width of the vehicle) 


Procedure: 


The function first computes useful computational parameters: 


P = 


Jjl 

Poo 


r 2 = (h-z cg )sin(X 2 +AX 2 r|) - (Li-x cg )cos(X 2 +Ax 2 T)) 


With these parameters equations (4.1-13) to (4.1-16), from [2], 
implemented to compute Xg, Zg, Mg, and Qe, respectively. 


Function Name: eomsphr 
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Description: 

This function computes the time derivatives of the twelve dynamic 



variables for time integration. 


Input Data: 



Symbol 
Text Code 

Description 

Units 

X 

xtot 

Body x-axis aeropropulsive force (per unit width of the vehicle) [lb/ft] 

Z 

ztot 

Body z-axis aeropropulsive force (per unit width of the vehicle) [lb/ft] 

M 

mtot 

Body y-axis aeropropulsive moment (per unit width of the 
vehicle) 

[ft-lb/ft] 

Qei 

qtot 

Generalized aeropropulsive force associated with the first 
generalized elastic degree of freedom (per unit width of the 
vehicle) 

[ft- lb/ft] 

h 

h 

Geometric altitude measured from the earth’s surface 

[ft] 

U 

u 

Vehicle forward speed 

[ft/s] 

V 

V 

Vehicle lateral speed 

[ft/s] 

w 

w 

Vehicle sink rate 

[ft/s] 

q 

q 

Vehicle rigid body pitch rate 

[rad/s] 

q 

eta 

First generalized elastic co-ordinate 

[-] 

q 

etadot 

First generalized elastic co-ordinate rate 

[-] 

X 

lam 

Latitude of vehicle position 

[rad] 

X 

tau 

Longitude of vehicle position 

[rad] 

Pi 

bl 

First quaternion component 

[-] 

P2 

b2 

Second quaternion component 

[-] 

P3 

b3 

Third quaternion component 

[-] 

P4 

b4 

Fourth quaternion component 

[-] 

<*>e 

ome 

Earth rotation rate 

[rad/s] 

Re 

re 

Earth radius 

[ft] 

go 

go 

Gravitational acceleration at Earth’s surface 

[ft/s 2 ] 

m 

m 

Vehicle mass (per unit width of the vehicle) 

[slug/ft] 

lyy 

lyy 

Vehicle moment of inertia about Body y axis (per unit width [ft 2 -slug/ft] 
of the vehicle) 

©i 

oml 

Undamped natural frequency of 1 st generalized elastic degree 
of freedom 

[rad/s] 

Cl 

zetl 

Damping ratio of 1 st generalized elastic degree of freedom 

[-] 

m El 

mel 

Generalized mass associated with the first generalized elastic [ft 2 -slug/ft] 
degree of freedom (per unit width of the vehicle) 
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Output Data: 


Symbol 
Text Code 

Description 

Units 

h 

dh 

Time rate of change of geometric altitude 

[ft/s] 

li 

du 

Body x-axis acceleration 

[ft/s 2 ] 

w 

dw 

Body z-axis acceleration 

[ft/s 2 ] 

q 

dq 

Body y-axis angular acceleration 

[rad/s 2 ] 

T] 

deta 

First generalized elastic co-ordinate rate 

[1/s] 

ii 

ddeta 

First generalized elastic co-ordinate acceleration 

[1/s 2 ] 

i 

dlam 

Time rate of change of vehicle latitude 

[rad/s] 

T 

dtau 

Time rate of change of vehicle longitude 

[rad/s] 

Pi 

dbl 

Time rate of change of first quaternion component 

[1/s] 

p2 

db2 

Time rate of change of second quaternion component 

[1/s] 

P3 

db3 

Time rate of change of third quaternion component 

[1/s] 

P4 

db4 

Time rate of change of fourth quaternion component 

[1/S] 

Procedure: 

The function begins by first computing the Vehicle-carrying frame 


to Body frame co-ordinate transformation matrix. Vehicle distance 
from the Earth's center is then computed. With these quantities 
equations (3.4) to (3.16) in Chapter 3 are implemented. The last 
four differential equations are implemented after the quantities p v , 
q v and r v are computed. Note: The twelve output quantities above 
are integrated and then fed into the "states" block (Figure Al) in 
the order that they appear above. Thus, for example, if the time 
history of w after a simulation is desired, it is found in the third 
column of the workspace variable "states." In addition, the 
quantities ti, w, q and i) are fed into the block "acc" in that order. 


B. PROCEDURAL OUTLINE OF ASUHS1 


The software development of ASUHS1 is centered around the "user-friendly" 
environment of SIMULINK®'s block diagram structure. It was designed to execute on 
any Apple Macintosh computer with SIMULINK® (or SIMULAB®) installed. . The 
suggested RAM size (1 Mbyte) allocated for MATLAB/S execution should be expanded 
to 2 Mbytes. This will allow the simulation to perform over 3,500 integration time steps 
under the current ASUHS1 set-up. The complete ASUHS1 package resides on a 3 1 /2" 
floppy disk which may be acquired from Professor David K. Schmidt at the Aerospace 
Research Center. 

The simulation package consists of seventeen MATLAB® functions and two script files. 
One of the functions is ASUHS1 itself, fifteen are functions called by ASUHS1 during 
execution and the other function is an auxiliary routine used to trim the vehicle (to be 
discussed later in this appendix). The two script files perform the task of loading vehicle 
parameters and initializing the controls and vehicle dynamic variables into the 
MATLAB® workspace. 

To initialize ASUHS1 execute the following steps: 

1) Load the script file inpar. 

2) Load the script file initstate . 

3) Activate the SIMULINK® block structure window and open ASUHS1. 

Each of these steps is detailed below. 

1) Within the script file inpar the Earth and atmosphere model parameters are specified, 
as well as vehicle and engine parameters. To load the parameter values presented in 
Section 6.1, simply type inpar from within MATLAB®'s command window: 

» i npar 

To change any of these parameter values the user may edit inpar, or create a new 
script file. To load the new parameter values into the workspace, follow the same 
procedure discussed above. (Note that if any parameters are altered, it may be 
necessary to determine a new trimmed flight condition.) 


2) ASUHS 1 initial conditions are specified within the script file initstate. This includes 
the vehicle dynamic variables and control inputs used to establish a trimmed flight 
condition. To load the initial conditions presented in Section 6.2, simply type 
initstate from within MATLAB®'s command window: 

ainitstate 

These initial conditions coincide with the parameter values in Section 6. 1 to provide a 
trimmed flight condition. As discussed in Section 6.2, this flight condition is 
obtained by specifying a vehicle Mach number, flight path angle, position and 
orientation while holding all angular rates and accelerations at zero. If a different 
flight condition is desired, the user must follow a trim procedure, such as the one 
outlined in Appendix C, to obtain a new trimmed flight condition. After obtaining 
new initial conditions, the file initstate may be edited or a new script file created to 
store the appropriate initial conditions. To load the new flight conditions into the 
workspace, follow the same procedure discussed above. 

Note that because some of the variables computed in initstate depend on values given by 
inpar, step (1) above should always precede step (2). 

3) To load ASUHS 1 in its block diagram structure, SIMULINK® s block diagram 
window must first be activated. To activate, type simulink within MATLAB®'s 
command window: 

»s i mu I ink 

SIMULINK®'s block library will appear. With the simulink window active (title bar 
in the simulink window is highlighted), open the file titled ASUHS 1. After a few 
seconds the ASUHS 1 window will appear as in Figure Al. 

The ASUHS 1 window displays the three control inputs, multiplexed together, as the input 
to the hypersonic vehicle model, block HV MODEL. The outputs of this block are the 
freestream and engine thermodynamic conditions, vehicle dynamic variables, trim 
accelerations, aeropropulsive forces, engine thrust and angle of attack. The last two 
outputs are provided for feedback control while the others are written to MATLAB®'s 
workspace for further analysis. Note that a time signal is also written to the workspace as 
a time reference for the other output data. 



Clock To Workspace 


Figure Al. ASUHS 1 block diagram in SIMULINK®. 

From any window within the ASUHS 1 block diagram structure the simulation parameters 
(start and stop time, integration algorithm, etc...) may be altered, via the "Parameters..." 
selection within the Simulation pull-down menu. The default algorithm used is RK-45 
with integration time step € [10' 5 , 10 -2 ] and relative error of 10' 3 . 

Given that the model parameters, initial conditions, ASUHS 1 and simulation parameters 
are loaded correctly, the simulation is now ready for execution. To run a simulation, 
select "Start" from the Simulation pull-down menu. The simulation will run up to the 
stop time. The responses in Figure 5, obtained for the numerical data presented in 
Chapter 6, are reproduced below in Figure A2. These plots may be used to verify proper 
execution of ASUHS 1. Due to the highly unstable behavior of this vehicle's 
configuration (lifting forebody), only short open-loop simulations are practical. 
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Figure A2. Simulation time responses to a doublet input in 5. 








C. TRIM PROCEDURES 


In this appendix a procedure is outlined which may be used to determine a trimmed flight 
condition for ASUHS 1. The procedure makes use of the nonlinear equation solver, /solve 
(included within MATLAB®), and relies on the flight condition structure given in Section 
6.2. If a different flight condition structure is desired, this procedure should be modified 
to account for the different fixed variable and free variable conditions. 

An equilibrium or trim condition for this flight condition structure requires that the 
accelerations u, w, q and rj be zero. The free parameters (a, 8, Ap and t|) are determined 
using the fsolve function, which is fully documented in [6], in conjunction with the coded 
function acczero. The latter function makes use of the aeropropulsive force and moment 
equations in the vehicle model, as well as the equations of motion. 

The computational trim procedure is as follows: For a given set of controls and dynamic 
variables, fsolve first employs acczero to compute the aeropropulsive forces and moments 
(via the equations in Chapter 4) and then the four accelerations of interest (via equations 
(3.4), (3.5), (3.7) and (3.8)). fsolve then compares each of the accelerations to zero and 
updates a, 8, Ad and T| accordingly. For each adjustment in the free parameters made by 
fsolve, acczero calculates the new accelerations. This iterative process continues until the 
four accelerations are zero. The inputs to fsolve are the function acczero itself and the 
first guesses for a, 8, Ap and T|. The outputs of fsolve are the values for a, 8, Ap and r\ 
which force the accelerations to zero. (These values are the trimmed vehicle and control 
settings and should be entered by the user into initstate for a 0 , 8^, Ap^ and T| 0 .) From 
MATLAB ®'s command window, type: 

acontro I s=fso I oe( ' acczero' , [-7 . 2 25.0 0.15 1.0]) 

where, -7.2°, 25.0°, 0.45 and 1.0 are the first guesses for a, 8, Ap and T|, respectively. 
(Note that all vehicle parameters and fixed variables are required to be specified within 
the acczero function. Thus any changes to the vehicle configuration or flight condition 
will require changes in the acczero code.) 

The critical question of how to choose the first guesses will now be addressed This is an 
important yet subtle step since using arbitrary values may produce unrealistic controls or 
vehicle configuration. A method used in this study plots the aeropropulsive forces and 
moments X, Z and M, as functions of angle of attack a for different pitch control surface 



settings 5 and different values of diffuser area ratio Ap. For each choice of a, 5 and Ap, 
the forces and moments are to be computed with rj trimmed (i.e., with rj=0). 

The method of trimming rj makes use of the spring restoring force in the vehicle elastic 
model. The approach begins by computing the aeropropulsive forces and moments for 
the given choice of a, 8 and Ap, and T|=l. The first generalized elastic co-ordinate 
acceleration Tj is then computed via equation (3.5) with T| set equal to zero. If rj is not 
zero, the elastic co-ordinate T| is perturbed by an amount At| so that the perturbed spring 
force, mEiajjAri (the model's torsional spring constant times the perturbed elastic co- 
ordinate), cancels this acceleration. The amount of perturbation is determined by 
examining equation (3.5) in Chapter 3, so that: 



Because the generalized force Qgi is dependent on the elastic deflection itself, several 
iterations (aeropropulsive force and acceleration computations followed by further 
perturbation of T|) are required to reduce this acceleration to zero. (Note that the final 
value of T| in this iterative procedure will likely be different for each choice of a, 8 and 
Ap. This final T| is not important; only the values of X, Z and M are of interest.) 

Along with the aeropropulsive forces and moments, the gravitational force components 
(vehicle weight) are also plotted. The control settings where the aeropropulsive force and 
gravitational force components intersect provide a good first guess for a, 8 and Ap. In 
Figure A3 we see that with Ap~0.45, 8~ 25° and a~-7.2°, the trim conditions for X, Z 
and M are all nearly met, so these are good first guesses for the first three inputs to the 
f, solve command. The first guess for T|, the fourth input, is not critical. A realistic value 
of 1.0 has been chosen here, which turns out to be close to the actual trim value given in 
Section 6.2. 

This appendix has briefly described a procedure to determine trim conditions. 
Alternatively, the trim command in SIMULINK® could be used to obtain trimmed 
vehicle and control settings. However, with the present ASUHS 1 set-up, experience has 
exposed fundamental difficulties associated with the implementation and convergence of 
the trim command. Nevertheless, there may exist some vehicle configurations for which 
the SIMULINK® trim command will execute appropriately. Thus the user who is 
familiar with the trim command may use it to determine a trim condition for his particular 
configuration if he so desires. 




Figure A3. Aeropropulsive forces plotted as functions of angle of attack. Different 
pitch control surface and diffuser area control settings were used. 



